We model a decision impact pathway is for school gardens as a general intervention for sustainable children’s food environments in urban Hanoi, Vietnam (Whitney et al. 2024).
Conceptual model of school gardens as an intervention. Should urban Hanoi school boards invest time and money in creating school gardens? Should they invest in formal STEM education as part of these gardens?
Simulation of the school garden intervention options:
# Source our model
source("CODAS_Garden_Model.R")
# Ensure consistent results with the random number generator
# not for each 'run' of the MC simulation but for
# consistency each time we run the entire simulation
set.seed(42)
garden_simulation_results <- mcSimulation(
estimate = estimate_read_csv("data/inputs_school_garden.csv"),
model_function = school_garden_function,
numberOfModelRuns = 1e4, #run 10,000 times
functionSyntax = "plainNames"
)
The Net Present Value (i.e. current value of the future benefits) of the garden decision options over 5 years of the intervention. For public and private schools the STEM costs are considered to be in the same garden space but with the additional costs and benefits of a full STEM education program. All options are compared to the same years of using the land for something that is not related to the garden, i.e. as a playground or for parking. Here we plot the distribution for the decision and frame the projected NPV.
For public schools
source("functions/plot_distributions.R")
plot_distributions(mcSimulation_object = garden_simulation_results,
vars = c("NPV_garden_public_school",
"NPV_garden_STEM_public_school"),
method = 'smooth_simple_overlay',
base_size = 7,
x_axis_name = "Comparative NPV outcomes")
For private schools
source("functions/plot_distributions.R")
plot_distributions(mcSimulation_object = garden_simulation_results,
vars = c("NPV_garden","NPV_garden_STEM"),
method = 'smooth_simple_overlay',
base_size = 7,
x_axis_name = "Comparative NPV outcomes")
The same results as boxplots
For public schools
source("functions/plot_distributions.R")
plot_distributions(mcSimulation_object = garden_simulation_results,
vars = c("NPV_garden_public_school", "NPV_garden_STEM_public_school"),
old_names = c("NPV_garden_public_school", "NPV_garden_STEM_public_school"),
new_names = c("NPV public school garden", "NPV public school garden with STEM"),
method = "boxplot_density",
base_size = 7,
x_axis_name = "Comparative NPV outcomes")
For private schools
source("functions/plot_distributions.R")
plot_distributions(mcSimulation_object = garden_simulation_results,
vars = c("NPV_garden","NPV_garden_STEM"),
old_names = c("NPV_garden","NPV_garden_STEM"),
new_names = c("NPV private school garden","NPV private school with STEM"),
method = "boxplot_density",
base_size = 7,
x_axis_name = "Comparative NPV outcomes")
summary(garden_simulation_results$y$NPV_garden)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1420.0 -171.5 382.3 692.8 1247.1 10218.6
Summary of the savings for the formal STEM education garden option
summary(garden_simulation_results$y$NPV_garden_STEM)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4017.3 -488.5 128.9 346.2 982.8 10175.3
Total expected costs for a school garden
summary(garden_simulation_results$y$total_costs)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 87.33 199.98 435.41 398.83 514.87 1474.13
Total expected costs for a school garden with STEM education
summary(garden_simulation_results$y$total_costs_STEM)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 143.1 357.2 839.6 929.8 1252.3 5011.9
First year expected costs for a school garden
summary(garden_simulation_results$y$Cashflow_garden1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -606.721 -95.383 8.306 65.859 174.797 1869.821
First year expected costs for a school garden with STEM education
summary(garden_simulation_results$y$Cashflow_garden_STEM1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -943.72 -237.06 -121.76 -77.60 44.09 1772.05
Here we assess value of information with the multi_EVPI
function.
# Subset the outputs from the mcSimulation function (y) by selecting the correct variables be sure to run the multi_EVPI only on the variables that we want. Find them with names(garden_simulation_results$y)
mcSimulation_table <- data.frame(garden_simulation_results$x,
garden_simulation_results$y[1:9])
Value of information for the garden option (no STEM)
source("functions/multi_EVPI.R")
# first_out_var is the first result variable in the table, "NPV_garden" in our case.
# names(garden_simulation_results$y)
# evpi <- multi_EVPI(mc = mcSimulation_table, first_out_var = "NPV_garden")
# save as a local .csv (takes ~ a day to run this)
# save(evpi,file="data/data_evpi.Rda")
load("data/data_evpi.Rda")
# open from saved file (last model run)
Generate the plots for EVPI PLS and cashflows for public and private schools
#Value of information the garden intervention decision
source("functions/plot_evpi.R")
plot_evpi_garden <- plot_evpi(evpi, decision_vars = "NPV_garden")
# Value of information for the garden option with formal STEM education.
# using the results of the same multi_EVPI
plot_evpi_STEM <- plot_evpi(evpi, decision_vars = "NPV_garden_STEM")
# Value of information for the public school garden option with no formal STEM education.
# using the results of the same multi_EVPI
plot_evpi_public <- plot_evpi(evpi, decision_vars = "NPV_garden_public_school")
# Value of information for the public school garden option with formal STEM education.
# using the results of the same multi_EVPI
plot_evpi_public_STEM <- plot_evpi(evpi, decision_vars = "NPV_garden_STEM_public_school")
# Cashflow of the garden option without formal STEM education
# This will be the cost for public and private schools over the intervention.
source("functions/plot_cashflow.R")
plot_cashflow_garden <- plot_cashflow(mcSimulation_object = garden_simulation_results,
cashflow_var_name = "Cashflow_garden")
# Cashflow of the garden option with formal STEM education
source("functions/plot_cashflow.R")
plot_cashflow_STEM <- plot_cashflow(mcSimulation_object = garden_simulation_results,
cashflow_var_name = "Cashflow_garden_STEM")
# Projection to Latent Structures (PLS)
# We use Projection to Latent Structures (PLS) model to get some sense of the correlation strength and direction for model variables and our outcome variables.
# For passive education garden option
source("functions/pls_model.R")
pls_result <- pls_model(object = garden_simulation_results,
resultName = names(garden_simulation_results$y)[1], # the "NPV_garden"
ncomp = 1)
# read in the common input table
input_table <- read.csv("data/inputs_school_garden.csv")
# source the plot function
source("functions/plot_pls.R")
plot_pls_garden <- plot_pls(pls_result,
input_table = input_table, threshold = 0.9)
#For school garden with formal STEM education
pls_result_STEM <- pls_model(object = garden_simulation_results,
resultName = names(garden_simulation_results$y)[2], # the "NPV_garden_STEM"
ncomp = 1)
plot_pls_STEM <- plot_pls(pls_result_STEM, input_table = input_table, threshold = 0.9)
NPV for public schools
#plot NPV as boxplot, pls, and evpi, evpi share names with pls
#each to a row
#shared x axis (values) and y axis (names)
dist_plot_NPV_garden_public <- plot_distributions(
mcSimulation_object = garden_simulation_results,
vars = "NPV_garden_public_school",
old_names = c("NPV_garden_public_school"),
new_names = c("NPV public school garden"),
method = "boxplot_density",
base_size = 7)
dist_plot_NPV_garden_STEM_public <- plot_distributions(
mcSimulation_object = garden_simulation_results,
vars = "NPV_garden_STEM_public_school",
old_names = c("NPV_garden_STEM_public_school"),
new_names = c("NPV public school garden with STEM"),
method = "boxplot_density",
base_size = 7)
# install.packages("devtools")
# devtools::install_github("thomasp85/patchwork")
library(patchwork)
dist_plot_NPV_garden_public / dist_plot_NPV_garden_STEM_public
## Warning: The following aesthetics were dropped during statistical transformation: x.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: x.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
PLS for public schools
plot_pls_garden / plot_pls_STEM
EVPI for publis schools
plot_evpi_STEM / plot_evpi_garden
NPV for private schools
#plot NPV as boxplot, pls, and evpi, evpi share names with pls
#each to a row
#shared x axis (values) and y axis (names)
dist_plot_NPV_garden <- plot_distributions(
mcSimulation_object = garden_simulation_results,
vars = "NPV_garden",
old_names = c("NPV_garden"),
new_names = c("NPV private school garden"),
method = "boxplot_density",
base_size = 7)
dist_plot_NPV_garden_STEM <- plot_distributions(
mcSimulation_object = garden_simulation_results,
vars = "NPV_garden_STEM",
old_names = c("NPV_garden_STEM"),
new_names = c("NPV private school with STEM"),
method = "boxplot_density",
base_size = 7)
# "NPV_garden_public_school", "NPV_garden_STEM_public_school"),
# install.packages("devtools")
# devtools::install_github("thomasp85/patchwork")
library(patchwork)
dist_plot_NPV_garden /
dist_plot_NPV_garden_STEM
## Warning: The following aesthetics were dropped during statistical transformation: x.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: x.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
PLS for private schools
plot_pls_garden / plot_evpi_garden
EVPI for private schools
plot_pls_STEM / plot_evpi_STEM
Cash flow plots of the garden option without formal STEM education. These are the expected costs for public and private schools over the intervention.
plot_cashflow_garden / plot_cashflow_STEM
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 15 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 15 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 5 rows containing missing values or values outside the scale range
## (`geom_line()`).
These figures display the Pareto-optimal solutions, representing the best trade-offs among the objectives of biodiversity, child health, and economic return. By focusing on these Pareto-optimal points, the analysis highlights solutions where improvements in one objective cannot be achieved without some compromise in at least one other.
source("functions/plot_pareto.R")
private_pareto <- pareto_front(
economic_return_garden = garden_simulation_results$y$NPV_garden,
child_health_garden = garden_simulation_results$y$health,
biodiversity_garden = garden_simulation_results$y$biodiversity,
economic_return_STEM = garden_simulation_results$y$NPV_garden_STEM,
child_health_STEM = garden_simulation_results$y$health_STEM,
biodiversity_STEM = garden_simulation_results$y$biodiversity,
plot_return = "scatter"
)
ggplotly(private_pareto)
knitr::include_graphics("figures/private_pareto_scatter.png")
knitr::include_graphics("figures/private_pareto_surface.png")
Trade offs for private school pareto
source("functions/pareto_front_posthoc.R")
private_pareto_posthoc <- pareto_front_posthoc(
economic_return_garden = garden_simulation_results$y$NPV_garden,
child_health_garden = garden_simulation_results$y$health,
biodiversity_garden = garden_simulation_results$y$biodiversity,
economic_return_STEM = garden_simulation_results$y$NPV_garden_STEM,
child_health_STEM = garden_simulation_results$y$health_STEM,
biodiversity_STEM = garden_simulation_results$y$biodiversity
)
## Number of Pareto-optimal points for STEM option: 47
## Number of Pareto-optimal points for Garden option: 66
##
## Summary of Pareto-optimal points for STEM option:
## economic_return biodiversity child_health
## Min. : -576.6 Min. : 718.5 Min. : 7.671
## 1st Qu.: 884.9 1st Qu.:1437.3 1st Qu.:24.519
## Median : 2178.2 Median :1867.4 Median :31.754
## Mean : 2883.6 Mean :2212.3 Mean :31.578
## 3rd Qu.: 4239.7 3rd Qu.:3000.3 3rd Qu.:38.132
## Max. :10218.6 Max. :5452.9 Max. :64.952
##
## Summary of Pareto-optimal points for Garden option:
## economic_return biodiversity child_health
## Min. :-2564.8 Min. : 402.6 Min. : 7.671
## 1st Qu.: 983.1 1st Qu.:1112.0 1st Qu.:18.787
## Median : 2264.2 Median :1533.7 Median :26.716
## Mean : 2536.1 Mean :1623.1 Mean :28.389
## 3rd Qu.: 4089.9 3rd Qu.:2086.7 3rd Qu.:36.457
## Max. :10175.3 Max. :3775.8 Max. :64.952
##
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
##
## smiths
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## Aggregation function missing: defaulting to length
private_pareto_posthoc
## $stem_summary
## economic_return biodiversity child_health
## Min. : -576.6 Min. : 718.5 Min. : 7.671
## 1st Qu.: 884.9 1st Qu.:1437.3 1st Qu.:24.519
## Median : 2178.2 Median :1867.4 Median :31.754
## Mean : 2883.6 Mean :2212.3 Mean :31.578
## 3rd Qu.: 4239.7 3rd Qu.:3000.3 3rd Qu.:38.132
## Max. :10218.6 Max. :5452.9 Max. :64.952
##
## $garden_summary
## economic_return biodiversity child_health
## Min. :-2564.8 Min. : 402.6 Min. : 7.671
## 1st Qu.: 983.1 1st Qu.:1112.0 1st Qu.:18.787
## Median : 2264.2 Median :1533.7 Median :26.716
## Mean : 2536.1 Mean :1623.1 Mean :28.389
## 3rd Qu.: 4089.9 3rd Qu.:2086.7 3rd Qu.:36.457
## Max. :10175.3 Max. :3775.8 Max. :64.952
##
## $trade_off_plot
source("functions/plot_pareto.R")
public_pareto <- pareto_front(
economic_return_garden = garden_simulation_results$y$NPV_garden_public_school,
child_health_garden = garden_simulation_results$y$health,
biodiversity_garden = garden_simulation_results$y$biodiversity,
economic_return_STEM = garden_simulation_results$y$NPV_garden_STEM_public_school,
child_health_STEM = garden_simulation_results$y$health_STEM,
biodiversity_STEM = garden_simulation_results$y$biodiversity,
plot_return = "scatter"
)
ggplotly(public_pareto)
knitr::include_graphics("figures/public_pareto_scatter.png")
knitr::include_graphics("figures/public_pareto_surface.png")
Trade offs for public school pareto
source("functions/pareto_front_posthoc.R")
public_pareto_posthoc <- pareto_front_posthoc(
economic_return_garden = garden_simulation_results$y$NPV_garden_public_school,
child_health_garden = garden_simulation_results$y$health,
biodiversity_garden = garden_simulation_results$y$biodiversity,
economic_return_STEM = garden_simulation_results$y$NPV_garden_STEM_public_school,
child_health_STEM = garden_simulation_results$y$health_STEM,
biodiversity_STEM = garden_simulation_results$y$biodiversity
)
## Number of Pareto-optimal points for STEM option: 56
## Number of Pareto-optimal points for Garden option: 71
##
## Summary of Pareto-optimal points for STEM option:
## economic_return biodiversity child_health
## Min. :-576.6 Min. : 0 Min. : 0.00
## 1st Qu.:1222.1 1st Qu.:1092 1st Qu.:21.40
## Median :2839.3 Median :1469 Median :27.83
## Mean :2978.7 Mean :1823 Mean :28.17
## 3rd Qu.:5080.4 3rd Qu.:2199 3rd Qu.:35.67
## Max. :7138.2 Max. :5453 Max. :64.95
##
## Summary of Pareto-optimal points for Garden option:
## economic_return biodiversity child_health
## Min. :-2564.85 Min. : 151.4 Min. : 7.671
## 1st Qu.: 66.99 1st Qu.:1003.4 1st Qu.:19.026
## Median : 2078.60 Median :1404.3 Median :24.732
## Mean : 2151.91 Mean :1576.8 Mean :27.016
## 3rd Qu.: 3952.36 3rd Qu.:2064.8 3rd Qu.:35.411
## Max. : 6583.17 Max. :3775.8 Max. :64.952
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## Aggregation function missing: defaulting to length
public_pareto_posthoc
## $stem_summary
## economic_return biodiversity child_health
## Min. :-576.6 Min. : 0 Min. : 0.00
## 1st Qu.:1222.1 1st Qu.:1092 1st Qu.:21.40
## Median :2839.3 Median :1469 Median :27.83
## Mean :2978.7 Mean :1823 Mean :28.17
## 3rd Qu.:5080.4 3rd Qu.:2199 3rd Qu.:35.67
## Max. :7138.2 Max. :5453 Max. :64.95
##
## $garden_summary
## economic_return biodiversity child_health
## Min. :-2564.85 Min. : 151.4 Min. : 7.671
## 1st Qu.: 66.99 1st Qu.:1003.4 1st Qu.:19.026
## Median : 2078.60 Median :1404.3 Median :24.732
## Mean : 2151.91 Mean :1576.8 Mean :27.016
## 3rd Qu.: 3952.36 3rd Qu.:2064.8 3rd Qu.:35.411
## Max. : 6583.17 Max. :3775.8 Max. :64.952
##
## $trade_off_plot
Here we provide a summary of the garden intervention options. We do
this with a summary table of the simulation results. We show the
percentage of missing values as well as the mean, median and standard
deviation (SD) for each output of our model simulations. We use the
gt_plt_summary() from {gtExtras} and with options from
{svglite}. The table shows the name, the plot overview as well as the
number of missing values, the mean, median and the standard deviation of
the distribution for all variables that were fed into the model from our
input table of uncertainty values.
# Subset the outputs from the mcSimulation function (y) to summarize only on the variables that we want.
# names(garden_simulation_results$x)
mcSimulation_table_x <- data.frame(garden_simulation_results$x[4:7,21:30]) #, ,32:41,43:70,73:76 also of possible interest
gtExtras::gt_plt_summary(mcSimulation_table_x)
| mcSimulation_table_x | ||||||
| 4 rows x 10 cols | ||||||
| Column | Plot Overview | Missing | Mean | Median | SD | |
|---|---|---|---|---|---|---|
| equipment_cost | 0.0% | 13.8 | 13.5 | 1.3 | ||
| construction_cost | 0.0% | 72.8 | 85.7 | 32.1 | ||
| garden_designing_costs | 0.0% | 11.7 | 11.9 | 1.3 | ||
| teacher_training_cost | 0.0% | 17.4 | 17.7 | 5.5 | ||
| school_board_planning | 0.0% | 7.1 | 7.2 | 1.1 | ||
| teaching_equipment | 0.0% | 9.6 | 9.5 | 1.9 | ||
| compost_starting | 0.0% | 8.3 | 8.4 | 1.9 | ||
| worm_starting | 0.0% | 6.5 | 5.7 | 2.2 | ||
| livestock_establishment_costs | 0.0% | 20.1 | 21.1 | 2.7 | ||
| fishpond_cost | 0.0% | 8.6 | 8.7 | 1.0 | ||
# a summary table with missing, mean, median and sd
The table shows the variable name, the plot overview as well as the number of missing values, the mean, median and the standard deviation of the distribution for variables that calculated in the model.
The full repository can be accessed at https://github.com/CWWhitney/urban_school_gardens